# ============================================================================
# CPU-GPU精度对比测试 - PlasticDamageConcrete3d (100单元)
# ============================================================================

puts "=========================================="
puts "CPU-GPU精度对比测试"
puts "网格规模: 10x10 = 100 Quad单元"
puts "材料: PlasticDamageConcrete3d (C30混凝土)"
puts "=========================================="

# ============================================================================
# 定义建模函数（复用）
# ============================================================================

proc buildModel {} {
    wipe
    model basic -ndm 2 -ndf 2

    # 网格参数
    set nx 80;  # X方向单元数
    set ny 80;  # Y方向单元数
    set Lx 5.0;  # X方向总长度 (m)
    set Ly 5.0;  # Y方向总长度 (m)

    set dx [expr $Lx / $nx];  # X方向单元尺寸
    set dy [expr $Ly / $ny];  # Y方向单元尺寸

    puts "网格尺寸: ${Lx}m x ${Ly}m (${nx}x${ny} = [expr $nx*$ny] 单元)"

    # C30混凝土材料
    nDMaterial PlasticDamageConcrete3d 1 30.0e9 0.2 3.0e6 30.0e6 0.6 0.5 2.0 0.75

    # 创建节点
    set nodeTag 1
    for {set j 0} {$j <= $ny} {incr j} {
        for {set i 0} {$i <= $nx} {incr i} {
            set x [expr $i * $dx]
            set y [expr $j * $dy]
            node $nodeTag $x $y
            incr nodeTag
        }
    }

    # 创建Quad单元
    set thickness 0.2
    set eleTag 1
    for {set j 0} {$j < $ny} {incr j} {
        for {set i 0} {$i < $nx} {incr i} {
            set n1 [expr $j * ($nx+1) + $i + 1]
            set n2 [expr $n1 + 1]
            set n3 [expr $n1 + ($nx+1) + 1]
            set n4 [expr $n1 + ($nx+1)]
            element quad $eleTag $n1 $n2 $n3 $n4 $thickness "PlaneStress" 1
            incr eleTag
        }
    }

    # 边界条件: 底部完全固定
    for {set i 1} {$i <= [expr $nx+1]} {incr i} {
        fix $i 1 1
    }

    # 左右侧边界: X方向固定
    for {set j 1} {$j < $ny} {incr j} {
        set leftNode [expr $j * ($nx+1) + 1]
        set rightNode [expr $leftNode + $nx]
        fix $leftNode 1 0
        fix $rightNode 1 0
    }

    # 顶部节点: X方向固定
    set topRowStart [expr $ny * ($nx+1) + 1]
    for {set i 0} {$i <= $nx} {incr i} {
        set nodeId [expr $topRowStart + $i]
        fix $nodeId 1 0
    }

    # 施加均布载荷（顶部节点）
    set totalLoad -50.0e3;  # 50 kN压缩载荷
    set nodeLoad [expr $totalLoad / ($nx + 1)]

    pattern Plain 1 Linear {
        for {set i 0} {$i <= $nx} {incr i} {
            set nodeId [expr $topRowStart + $i]
            load $nodeId 0.0 $nodeLoad
        }
    }

    puts "模型构建完成"
    puts "  节点数: [expr ($nx+1) * ($ny+1)]"
    puts "  单元数: [expr $nx * $ny]"
    puts "  顶部总载荷: [expr -$totalLoad/1000] kN"

    # 返回关键参数
    return [list $nx $ny $topRowStart]
}

# ============================================================================
# 定义监测节点选择函数
# ============================================================================

proc selectMonitorNodes {nx ny topRowStart} {
    # 选择5个关键节点进行监测
    # 1. 左下角 (固定节点，应该位移为0)
    # 2. 右下角 (固定节点)
    # 3. 左上角
    # 4. 右上角
    # 5. 中心顶部

    set bottomLeft 1
    set bottomRight [expr $nx + 1]
    set topLeft $topRowStart
    set topRight [expr $topRowStart + $nx]
    set topCenter [expr $topRowStart + $nx/2]

    return [list $bottomLeft $bottomRight $topLeft $topRight $topCenter]
}

# ============================================================================
# CPU分析
# ============================================================================

puts "\n=========================================="
puts "第一阶段: CPU分析"
puts "=========================================="

# 构建模型
set params [buildModel]
set nx [lindex $params 0]
set ny [lindex $params 1]
set topRowStart [lindex $params 2]

# 选择监测节点
set monitorNodes [selectMonitorNodes $nx $ny $topRowStart]

puts "\n配置CPU求解器..."
constraints Plain
numberer Plain
system BandGeneral
test NormDispIncr 1.0e-6 20
algorithm Newton
integrator LoadControl 0.1

analysis Static

puts "开始CPU分析 (10步)..."
set startTime [clock milliseconds]

set numSteps 10
set ok 0

for {set i 1} {$i <= $numSteps} {incr i} {
    set ok [analyze 1]
    if {$ok != 0} {
        puts "❌ CPU分析在第${i}步失败"
        break
    }
    if {$i % 5 == 0} {
        puts "  Step $i/$numSteps 完成"
    }
}

set endTime [clock milliseconds]
set cpuTime [expr $endTime - $startTime]

if {$ok == 0} {
    puts "✓ CPU分析完成，耗时: $cpuTime ms"

    # 保存CPU结果
    set outFile [open "cpu_results.txt" w]
    puts $outFile "=== CPU分析结果 ==="
    puts $outFile "网格: ${nx}x${ny} = [expr $nx*$ny] 单元"
    puts $outFile "分析步数: $numSteps"
    puts $outFile "总时间: $cpuTime ms"
    puts $outFile "平均时间/步: [format "%.2f" [expr $cpuTime / double($numSteps)]] ms"
    puts $outFile ""
    puts $outFile [format "%-10s %-15s %-15s" "Node" "Disp_X(m)" "Disp_Y(m)"]
    puts $outFile "--------------------------------------------"

    # 保存CPU位移结果到数组
    global cpuDisp
    foreach nodeId $monitorNodes {
        set ux [nodeDisp $nodeId 1]
        set uy [nodeDisp $nodeId 2]
        set cpuDisp($nodeId,x) $ux
        set cpuDisp($nodeId,y) $uy
        puts $outFile [format "%-10d %-15.6e %-15.6e" $nodeId $ux $uy]
    }

    close $outFile
    puts "CPU结果已保存到: cpu_results.txt"
} else {
    puts "❌ CPU分析失败，无法继续对比测试"
    exit
}

# ============================================================================
# GPU分析
# ============================================================================

puts "\n=========================================="
puts "第二阶段: GPU分析"
puts "=========================================="

# 重新构建模型
set params [buildModel]
set nx [lindex $params 0]
set ny [lindex $params 1]
set topRowStart [lindex $params 2]

puts "\n配置GPU求解器..."
constraints Plain
numberer Plain
system GPUSolver -debugLevel 1 -enableModule all
test NormDispIncr 1.0e-6 20
algorithm Newton
integrator LoadControl 0.1

analysis Static

puts "开始GPU分析 (10步)..."
set startTime [clock milliseconds]

set ok 0

for {set i 1} {$i <= $numSteps} {incr i} {
    set ok [analyze 1]
    if {$ok != 0} {
        puts "❌ GPU分析在第${i}步失败"
        break
    }
    if {$i % 5 == 0} {
        puts "  Step $i/$numSteps 完成"
    }
}

set endTime [clock milliseconds]
set gpuTime [expr $endTime - $startTime]

if {$ok == 0} {
    puts "✓ GPU分析完成，耗时: $gpuTime ms"

    # 保存GPU结果
    set outFile [open "gpu_results.txt" w]
    puts $outFile "=== GPU分析结果 ==="
    puts $outFile "网格: ${nx}x${ny} = [expr $nx*$ny] 单元"
    puts $outFile "分析步数: $numSteps"
    puts $outFile "总时间: $gpuTime ms"
    puts $outFile "平均时间/步: [format "%.2f" [expr $gpuTime / double($numSteps)]] ms"
    puts $outFile ""
    puts $outFile [format "%-10s %-15s %-15s" "Node" "Disp_X(m)" "Disp_Y(m)"]
    puts $outFile "--------------------------------------------"

    # 保存GPU位移结果到数组
    global gpuDisp
    foreach nodeId $monitorNodes {
        set ux [nodeDisp $nodeId 1]
        set uy [nodeDisp $nodeId 2]
        set gpuDisp($nodeId,x) $ux
        set gpuDisp($nodeId,y) $uy
        puts $outFile [format "%-10d %-15.6e %-15.6e" $nodeId $ux $uy]
    }

    close $outFile
    puts "GPU结果已保存到: gpu_results.txt"
} else {
    puts "❌ GPU分析失败，无法继续对比测试"
    exit
}

# ============================================================================
# 精度对比分析
# ============================================================================

puts "\n=========================================="
puts "第三阶段: 精度对比分析"
puts "=========================================="

global cpuDisp gpuDisp

set outFile [open "cpu_gpu_comparison.txt" w]
puts $outFile "=== CPU-GPU精度对比报告 ==="
puts $outFile "网格规模: ${nx}x${ny} = [expr $nx*$ny] 单元"
puts $outFile "材料: PlasticDamageConcrete3d (C30混凝土)"
puts $outFile ""
puts $outFile "性能对比:"
puts $outFile "  CPU总时间: $cpuTime ms (平均 [format "%.2f" [expr $cpuTime / double($numSteps)]] ms/步)"
puts $outFile "  GPU总时间: $gpuTime ms (平均 [format "%.2f" [expr $gpuTime / double($numSteps)]] ms/步)"
set speedup [expr double($cpuTime) / double($gpuTime)]
puts $outFile "  加速比: [format "%.2f" $speedup]x"
puts $outFile ""
puts $outFile "位移精度对比:"
puts $outFile [format "%-10s %-15s %-15s %-15s %-15s %-12s %-12s" "Node" "CPU_Ux(m)" "CPU_Uy(m)" "GPU_Ux(m)" "GPU_Uy(m)" "Error_x(%)" "Error_y(%)"]
puts $outFile "----------------------------------------------------------------------------------------------------------------------------"

# 打印表头到控制台
puts ""
puts [format "%-10s %-15s %-15s %-15s %-15s %-12s %-12s" "Node" "CPU_Ux(m)" "CPU_Uy(m)" "GPU_Ux(m)" "GPU_Uy(m)" "Error_x(%)" "Error_y(%)"]
puts "----------------------------------------------------------------------------------------------------------------------------"

# 计算误差统计
set maxErrorX 0.0
set maxErrorY 0.0
set sumErrorX 0.0
set sumErrorY 0.0
set count 0

foreach nodeId $monitorNodes {
    set cpuUx $cpuDisp($nodeId,x)
    set cpuUy $cpuDisp($nodeId,y)
    set gpuUx $gpuDisp($nodeId,x)
    set gpuUy $gpuDisp($nodeId,y)

    # 计算相对误差 (%)
    # 注意：如果CPU值接近0，使用绝对误差
    if {abs($cpuUx) > 1.0e-15} {
        set errorX [expr abs(($gpuUx - $cpuUx) / $cpuUx) * 100.0]
    } else {
        set errorX [expr abs($gpuUx - $cpuUx) * 100.0]
    }

    if {abs($cpuUy) > 1.0e-15} {
        set errorY [expr abs(($gpuUy - $cpuUy) / $cpuUy) * 100.0]
    } else {
        set errorY [expr abs($gpuUy - $cpuUy) * 100.0]
    }

    # 更新最大误差
    if {$errorX > $maxErrorX} {set maxErrorX $errorX}
    if {$errorY > $maxErrorY} {set maxErrorY $errorY}

    # 累加误差
    set sumErrorX [expr $sumErrorX + $errorX]
    set sumErrorY [expr $sumErrorY + $errorY]
    incr count

    # 输出到文件和控制台
    set line [format "%-10d %-15.6e %-15.6e %-15.6e %-15.6e %-12.6f %-12.6f" $nodeId $cpuUx $cpuUy $gpuUx $gpuUy $errorX $errorY]
    puts $outFile $line
    puts $line
}

# 计算平均误差
set avgErrorX [expr $sumErrorX / double($count)]
set avgErrorY [expr $sumErrorY / double($count)]

puts ""
puts $outFile ""
puts $outFile "误差统计:"
puts $outFile "  最大相对误差 (X方向): [format "%.6e" $maxErrorX]%"
puts $outFile "  最大相对误差 (Y方向): [format "%.6e" $maxErrorY]%"
puts $outFile "  平均相对误差 (X方向): [format "%.6e" $avgErrorX]%"
puts $outFile "  平均相对误差 (Y方向): [format "%.6e" $avgErrorY]%"

set maxError [expr $maxErrorX > $maxErrorY ? $maxErrorX : $maxErrorY]
set threshold 1.0e-4;  # 0.0001% = 1e-6相对误差

if {$maxError < $threshold} {
    puts $outFile "  ✓ 精度验证通过 (最大误差 < [format "%.6e" $threshold]%)"
    set status "通过"
} else {
    puts $outFile "  ⚠ 精度误差较大 (最大误差 > [format "%.6e" $threshold]%)"
    set status "警告"
}

close $outFile

puts "=========================================="
puts "精度对比统计:"
puts "  最大相对误差 (X): [format "%.6e" $maxErrorX]%"
puts "  最大相对误差 (Y): [format "%.6e" $maxErrorY]%"
puts "  平均相对误差 (X): [format "%.6e" $avgErrorX]%"
puts "  平均相对误差 (Y): [format "%.6e" $avgErrorY]%"
puts "  验证状态: $status"
puts "=========================================="

puts "\n性能对比:"
puts "  CPU时间: $cpuTime ms"
puts "  GPU时间: $gpuTime ms"
puts "  加速比: [format "%.2f" $speedup]x"
puts "=========================================="

puts "\n完整对比报告已保存到: cpu_gpu_comparison.txt"

wipe
